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ABSTRACT 

We present an initial report on 160 simulations of a highly simplified model of 
the post-bounce core-collapse supernova environment in three spatial dimensions 
(3D). We set different values of a parameter characterizing the impact of nuclear 
dissociation at the stalled shock in order to regulate the post-shock fluid velocity, 
thereby determining the relative importance of convection and the stationary 
accretion shock instability (SASI). While our convection-dominated runs comport 
with the paradigmatic notion of a ‘critical neutrino luminosity’ for explosion at 
a given mass accretion rate (albeit with a nontrivial spread in explosion times 
just above threshold), the outcomes of our SASI-dominated runs are much more 
stochastic: a sharp threshold critical luminosity is ‘smeared out’ into a rising 
probability of explosion over a ~ 20% range of luminosity. We also find that 
the SASI-dominated models are able to explode with 3 to 4 times less efficient 
neutrino heating, indicating that progenitor properties, and fluid and neutrino 
microphysics, conducive to the SASI would make the neutrino-driven explosion 
mechanism more robust. 

Subject headings: hydrodynamics — instabilities — stars: evolution — super¬ 
novae: general 
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1. Introduction 


Despite a notorious lack of robustness, the delayed neutrino-driven mechanism has per¬ 
sisted for 30 years as the leading paradigm of typical core-collapse supernova explosions 
( Bethe h Wilson||1985| Kotake et al. 2012; Janka 2012; Burrows 2013). Gravitational col¬ 
lapse of a massive (> 8 Mq) star halts with a ‘bounce’ at nuclear density, launching a 

shock wave that stalls at 100 — 200 km, a few times the radius of the hot nascent neutron 
star. Neutrino fluxes from the neutron star heat the outer post-shock region. In spherically 
symmetric distillations of the problem, if the neutrino luminosity is sufficiently high for a 
given mass accretion rate, the ram pressure of infalling matter is overcome and the shock re¬ 
sumes its outward journey (Burrows & Goshy 1993; Janka||2001i Yamasaki &: Yamada||2007 


Pejcha fc Thompson||2012 i Fernandez||2012 ; Muller fc Janka||2015 ; Murphy &: DolerL^|2015 ). 
Multidimensional fluid instabilities increase the dwell time of matter in the heating region 
and the average stalled shock radius relative to spherical symmetry; these effects are often 


characterized as 

reducing this critical or threshold luminosity for explosion (Burrows et al. 

1995; 

Janka & Mueller] 1996; 

Yamasaki & Yamada||2005; lOhnishi et aL|2006; Yamasaki & 

Yamada||2006 ; Iwakami et al. 

2008; [Murphy & Burrows||2008 ; Fernandez & Thompson||2009; 

Hanke et al. 2012; Dolence et al. |2013[ Gouch|2013; Gouch & Ott 2015; IMiiller & Janka 

2015; 

Fernandez 

2015). Beyond these generally more favorable conditions, it seems that the 


runaway expansion of one or a few large, persistent, buoyant bubbles, inflated by neutrino 


heating, is the final proximate agent propelling the shock on its way (Herant 1995; Thomp- 


son 2000| [Fernandez & Thompson 2009| [Dolence et al. 2013| [Wongwathanarat et al.([2013| 

Gouch|2013; [Handy et al.|2014l [Fernandez et al.[[2014 

; Muller & Janka|2015l Fernandez|2015; 

Lentz et al. 

2015 

). Here we consider two questions related to the delayed neutrino-driven 


mechanism. 

First, does this system exhibit qualitatively significant sensitive dependence on initial 
conditions? How reproducible and representative is a single simulation, compared to others 
with nearly identical initial states? The concern is that 3D simulations with more sophisti¬ 
cated neutrino transport take tens of millions of processor hours and months to run; only a 
handful have been reported (Hanke et al.||2013; Tamborra et al. 2014; Melson et al.||2015b| 


Lentz et al. 2015). 


Second, might the stationary accretion instability (SASI) materially contribute to su¬ 
pernova explosions (Hanke et al.||2012; Muller et al.||2012| Hanke et al.||2013; Fernandez et al. 


; Tamborra et al.[[2014; 

Iwakami et al.[[2014; 

Melson et al. [[2015a; 

Fernandez 

2015 

); or is 


it generically subdominant to an extent that convection is all that matters ( 

Burrows et al. 

20121 Murphy et aL||2013; Dolence et aL||2013HHandy et aL||2014; ' 

Fakiwaki et aL||2014; Couch 

& O’Gonnor 2014; Abdikamalov et al. 2015; Melson et al. 2015b 

; Lentz et al. 

2015)? In the 
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absence of heating, a stationary accretion shock does not remain stationary: small aspherici- 
ties result in large-scale (multipole £= 1,2), nonlinear, oscillatory motions—typically, spiral 
in 3D ( Blondin et al.|[^03| [Blondin fc Mezzacappa 2007; Fernandez 2010; Guilet & 


waves 


Foghzzo||2012). This increases the average shock radius, but not to explosion; neutrino heat¬ 


ing, which drives convection and inflates large buoyant bubbles, is ultimately required. But 
strong neutrino heating tends to suppress the SASI, as its persistent larger-scale bubbles 


interrupt cyclic, globally coherent SASI flows (Burrows et al. 2012; Fernandez et al. 2014 


Couch & O’Connor 2014; Iwakami et al. 2014; Fernandez 2015). The question is whether 


the SASI facilitates the heating and inflation of large buoyant bubbles to any important ex¬ 


tent, by lengthening the advection timescale and generating entropy gradients (Scheck et al. 


2008). 


2. Model 


We explore these questions with a simplified model of the post-bounce supernova envi¬ 
ronment, closely following Fernandez et al. (2014) and Fernandez (2015), allowing for many 
realizations and control over the nature of the dominant instability 


We obtain initial conditions from stationary, spherically symmetric, and nonrelativistic 
fluid equations for baryon, momentum, and energy conservation. We use an ideal gas equa¬ 
tion of state with adiabatic index 7 = 4/3, central gravity with fixed point mass M = 1.3 Mq, 
and constant mass accretion rate M = 0.3 Mq A parameter £ mimicks nuclear dissoci¬ 
ation in the shock jump conditions. We use the energy source 


g,= (^-Ap^/^jpe(Mo-M), 


( 1 ) 


where B and A parametrize the magnitude of neutrino heating and cooling, r is the radius, 
p is the pressure, p is the mass density, A4 is the Mach number, and A4o = 2 is a reference 
value separating the pre- and post-shock regions; the Heaviside step function 0 restricts 
heating and cooling to the post-shock region. 


Initial conditions in the pre- and post-shock regions are joined by the shock jump con¬ 
ditions. An inner reflecting boundary at tns = 40 km represents the (fixed) surface of the 
hot nascent neutron star. We choose a shock radius rso = 100 km in the absence of neutrino 
heating (B = 0). A vanishing Bernoulli parameter and A4 = 6 at rgo determine the nearly 
free-fall pre-shock solution. The e-dependent shock jump conditions determine the initial 
conditions of integration from the shock to the inner boundary. All simulations with a given 
£ use the same A. We determine these A values empirically, as those which yield zero velocity 
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at tns with B = 0. In setting up initial profiles for increasing values of B (with e and A 
fixed), the shock radius rs > rso is empirically determined as that which again yields zero 
velocity at tns- 


We evolve the system with GenASiS (Cardall et al. 2014; Cardall &: Budiardja||2015] ), 
using a multilevel mesh in Cartesian coordinates. The coarsest level covers [—640 km, +640 km] 
in 3D. Three additional levels consist of nested spheres of radius 320 km, 160 km, and 80 km. 
The coarsest contains 128^ cells (cell width 10 km); successive levels are refined by a factor of 
2 (finest level cell width 1.25 km). We interpolate the spherically symmetric initial conditions 
to the 3D multilevel mesh and add random 0.1% pressure perturbations. We implement the 
mock-up of nuclear dissociation with an additional continuity equation for the density of 
‘dissociated baryons,’ with appropriate sources for this and the energy equation. 


3. Results 


We present ‘convection-dominated’ (C-series) and ‘SASI-dominated’ (S-series) runs. 
Convection only develops from small perturbations if the post-shock infall velocity is suf¬ 
ficiently slow; otherwise, incipient buoyant perturbations are advected to the neutron star 
before they can appreciably grow. In that case, the SASI is more free to operate. The 
Foglizzo number y—the ratio of advection time to buoyancy timescale in the region of net 
heating beneath the shock—is a quantitative measure of this, with y > 3 allowing convection 
to develop from linear perturbations (Foglizzo et al. 2006). 


Our C-series and S-series differ in dissociation parameter e, which determines their 
convection-dominated vs. SASI-dominated character. By removing energy from the fluid, 
larger values increase the shock compression ratio and reduce the post-shock infall speed; 
this increases the advection time scale, and therefore y. The values of y in Table for 
the initial profiles used in the C-series and S-series are respectively well above and below 
the value y == 3 separating the convection-dominated and SASI-dominated regimes. The 
perhaps extreme £ values serve to separate these regimes; the messier truth may he between. 


We summarize our simulations in Table and more selectively in Figures and For 
each series, we run 10 simulations for each value of B in Table [T| which summarizes some 
parameters of the initial conditions and statistics of the outcomes. Figures [T] and survey 
the differing shock behavior of each series as a function of heating parameter: for selected 
values of B below, at, and above ‘threshold’, the £ = 0,1 real spherical harmonic coefficients 
ttirn of the shock surface are displayed for several simulations, normalized such that aoo is 
the average shock radius, and an, ai_i, and aio are the average x, y, and .2 positions of the 
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Table 1: Simulation summary. 


Dominant 

6: 

A 

B 

rs 

X 

A400 

^400,min 

At 400 

Instability 

(^0) 

(A) 

(Bo) 

(rso) 



(s) 

(s) 

Convection 

0.3 

0.0131 

0.700 

1.202 

6.27 

0 



(C-Series) 



0.725 

1.212 

6.62 

0 






0.750 

1.223 

6.97 

0 






0.775 

1.234 

7.33 

0 






0.800 

1.245 

7.71 

10 

0.333 

0.180 




0.825 

1.257 

8.09 

10 

0.209 

0.076 




0.850 

1.269 

8.49 

10 

0.157 

0.019 




0.875 

1.282 

8.89 

10 

0.141 

0.010 

SASI 

0.0 

0.1419 

0.950 

1.188 

0.93 

0 



(S-Series) 



0.975 

1.195 

0.97 

1 

0.268 





1.000 

1.202 

1.02 

4 

0.241 





1.025 

1.210 

1.06 

4 

0.219 





1.050 

1.217 

1.11 

9 

0.206 





1.075 

1.225 

1.15 

9 

0.239 





1.100 

1.233 

1.20 

9 

0.243 





1.125 

1.241 

1.25 

10 

0.213 

0.603 


Note. — Parameters of the initial conditions ( 5 , A, and y), and the number (A^ 4 oo): earliest 

time (t 4 oo,min )5 and time spread (At 4 oo) of explosions obtained from 10 runs per value of heating pa¬ 
rameter B (where ‘explosion’ here means the average shock radius reaching 400 km). Reference values: 
So = GM/rso, where G is Newton’s constant and rso = 100 km; Aq = 145 MeV s“^, 

where rriu is the atomic mass unit, po = (lla^/12)(2 MeV)^, and = 8.56 x 10^^ MeV~^cm“^; and 
Bq = 160 my^(10^ km)^ MeV s“^ ( Janka|[200l| ). 
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Fig. 1.— C-Series shock radius decomposition coefficients for selected B values. Solid lines 
are duds; dot-dashed lines are explosions. Upper panels: average shock radius aoo for all 10 
runs for each value of B. Lower panels: average x, y, and z shock positions (on, ai_i, oio), 
with only half the runs and the first 0.5 s for clarity (the time scales differ from the upper 
panels). 



Fig. 2.— As in Figure [T} for the S-series. 
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shock (Burrows et al. 2012). 


For the C-series, the number of explosions A^4oo as a function of B tabulated in Table [T] 
comports with the paradigmatic notion of a ‘critical’ luminosity: A/400 jumps from 0 to 10 
between values of B separated by only ~ 3%. Mild stochasticity is evident in the spread of 
explosion times (At)4oo in the C-series runs, which decreases with increasing 5, as does the 
earliest explosion time ^ 400 ,min- 


However, the critical luminosity concept is challenged by the more extreme stochasticity 
of the S-series, which extends beyond explosions times to whether an explosion occurs at 
all (at least within our Is simulation tim^: this ‘threshold’ is smeared out, with N^qq in 
Table [T] increasing from 0 to 10 over a ~ 20% range of B. The spread of explosion times is 
only given for the largest B value, the only one for which all 10 runs explode. This spread 
(At)4oo = 0.603 s is signihcantly larger than in the C-series. Unlike the tabulated values 
of ^400,min for the C-series, there is no clear trend of earliest explosion time with increasing 
B] apparently, if an explosion can happen at all, it has a decent chance of occurring in the 
initial expansion, at roughly the same time for all tabulated B. 

Figures and [^reinforce the outcomes in Table Whereas the C-series runs in Figure [T] 
show either steady shock expansion (explosions) or a stationary radius (duds), the strong 
influence of the SASI yields large, cyclic excursions on roughly the SASI growth timescale 
in the upper panels of Figure]^ Again, the explosion threshold is not sharp for the S-series: 
the upper middle panel of Figure shows some runs exploding, and others not, within the 
time frame of the simulation. Note also a telltale signature of the SASI in the lower panels of 
Figure [2] —the oscillatory dipole coefficients. These grow more irregular and episodic when 
convection develops and interrupts their global coherence. Such oscillations are absent from 
the lower panels of Figure [Tj dipole coefficients in these convection-dominated models grow in 
a non-oscillatory manner as explosions take off and bubbles merge into dominant directions 
of expansion. 

Figure shows sample C-series and S-series runs. In the exploding convective model 
(left panels), Rayleigh-Taylor plumes, driven by a negative entropy gradient generated by 
neutrino heating, grow from small scales and merge into larger bubbles, which become suf¬ 
ficiently buoyant and persistent to expand the shock to our outer boundary. In the SASI- 
dominated example (right-panels), a hint of spiral behavior manifests in the upper left cross 


^With our stationary luminosity and accretion rate, the concept of a critical luminosity might be recovered 
if the simulations ran indefinitely. But in nature, luminosity declines with time, leaving a finite time for 
explosion to occur. Thus, if SASI-dominated explosions occur in nature, the stochasticity of outcomes 
reported here may plausibly hold. 







Fig. 3.— Left panels: Sample convection-dominated simulation with B = 0.800 5o at t = 
0.1s, and at t = 0.25 s as the explosion takes off. Right panels: Sample SASI-dominated 
simulation with B — 1.025 Rq at t = 0.2 s, and Ski t — 0.95 s as the explosion takes off. 
Surface plots of two entropy values show the shock (cyan) and some structure in the post¬ 
shock region (purple). Entropy cross sections through the origin are projected onto the back 
walls. The average, maximum, and minimum shock radii are shown in the lower left of each 
panel, with the solid and dotted portions of these traces indicating times before and after 
the moment shown. For each simulation, the zero of entropy is taken to be the immediate 
post-shock value at t = 0. 
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section of the upper right panel, with a plunging stream of low-density material around 11 
o’clock separating counterrotating flows; but the behavior is not pure SASI, as smaller-scale 
Rayleigh-Taylor plumes are also present. 


Animations of S-series runs suggest that the long-time-scale, quasi-periodic excursions 
of the shock radius result from the palpable influence of both the SASI and neutrino heating. 
At small shock radii, SASI spiral waves appear. These increase the shock radius and gen¬ 
erate entropy gradients to a point favorable to convection. Even without neutrino heating, 
the Rayleigh-Taylor instability is a SASI saturation mechanism (Guilet et al. 2010); here, 
neutrino heating takes buoyancy beyond saturation of the SASI to interruption of the SASI. 
When convective bubbles are sufficiently persistent to interrupt the global SASI oscillation, 
but not sufficiently buoyant to drive runaway shock expansion, the shock radius collapses 
and the cycle begins anew. In a final cycle towards the end of the run in the right panels 
of Figure!^ the actions of the SASI and convection coincide and reinforce each other to the 
point that a large, buoyant bubble drives runaway expansion. 


This interpretation is supported by the time series of SASTinduced angular momentum 
L, and y, in the upper and middle panels of Figure The radius changes in the lower left 
corners of the right panels of Figure]^ are manifestly reflected in |L| and y. Strikingly, the 
y peaks are close to the critical value y ~ 3 for convective activity to persist. The shock 
expansion occasioned by the SASI gives convection a chance to take hold; but if runaway 
expansion does not immediately result, the shock recedes to a point that the spiral SASI is 
renewed with a stochastic reorientation, as indicated by the L component traces. 


Finally, we compare the heating efficiencies in C-series and S-series runs near threshold, 
which suggests that SASTdominated systems are significantly ‘easier’ to explode. Computed 
from spherically averaged radial profiles, the heating efficiency rj is the integrated net heating 
rate in the ‘gain region’ (where the net heating is positive) divided by the sum of the 
luminosities of the neutrino species responsible for heating -|- which is related to 
our heating parameter B ( Janka|200T )]. The lower panel of Figure]^ shows rj for the C-series 
runs just above threshold {B = 0.800 Bq) and for the S-series runs towards the middle of 
the ‘smeared out’ threshold region [B = 1.025 Rq)- Only 4 of these SASTdominated runs 
explode, but the heating efficiencies of these exploding runs are not drastically different from 
the 6 duds. For explosions, we see that heating efficiencies 3 to 4 times higher are needed 
for convection-dominated runs than for SASTdominated runs. 


The lower efficiencies needed for SASTdominated systems to explode imply that they 
are more ‘forgiving,’ in that they require less from the neutrino heating. This is physically 
plausible. In convection-dominated systems, the neutrino heating does double duty: in addi¬ 
tion to inflating bubbles to runaway buoyancy, neutrino-driven convection must be initially 
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Fig. 4.— Upper and middle panels: Angnlar momentnm magnitude (solid) and components 
(dot-dashed), and Foglizzo number x, of the SASI-dominated run shown in the right panels of 
Figure]^ Lower panel: Heating efficiency rj of convection-dominated runs with B = 0.800 Bq 
(upper curves) and SASI-dominated runs with B = 1.025 5o (lower curves). 
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strong enough to expand the shock against gravity and ram pressure to an extent that makes 
such inflation possible; whereas in SASI-dominated systems, the SASI initially does much of 
the ‘heavy lifting’ in terms of initial shock expansion and generation of entropy gradients, 
such that all that is required of the neutrino heating is the inflation of bubbles in conditions 
that are already more favorable. 


4. Conclusion 


We ran 160 simulations with a highly simplified model of the post-bounce supernova 
environment in 3D. This simplified model allows control over the nature of the dominant 
multidimensional instability affecting shock evolution. 


Do core-collapse supernovae exhibit qualitatively significant sensitive dependence on ini- 


tial conditions? 

Wongwathanarat et al. 

(2013)^ 

, Handy et al. 

(2014), 

, Takiwaki et al. 

(2014), 

and 

Iwakami et al. 

(2014 

)—with more physics, though in most cases still parametrized, and 


with relative handfuls of 3D realizations—And quantitative spreads in explosion time and en¬ 
ergy, along with larger variations in ejecta morphology, neutron star spins and kick velocities, 
and distributions of nucleosynthesis products. We also see quantitative spreads in explosion 
time in our convection-dominated models, which decrease with increasing luminosity. 


However, we find that stochasticity extends to the qualitative outcome of explosion vs. 
dud in the SASI-dominated case. Out of 10 simulations for each value of heating parameter 
B listed for the SASI-dominated models in Table the number of explosions increases from 
0 to 10 over a ~ 20% range in luminosity. This stochasticity in qualitative outcome arises 
from the strong interplay of two instabilities. The SASI is enabled first, increasing the 
shock radius to a point that allows convection to begin. Subsequent evolution depends on 
whether neutrino heating drives runaway bubble inflation at this key moment. If not, the 
shock recedes and the cycle repeats with the SASI again taking over. This long-period cyclic 


interaction has been glimpsed before—see 

Iwakami et al. 

(2014 

) Figure 2, 

Fernandez 

(2015) 

Figure 3, and some self-consistent models ( 

Hanke et al. 

2013; 

Tamborra et al. 

2014 

). But 


a larger number of realizations and finer sampling of heating strength have here revealed 
that this cyclic interaction leads to more qualitative stochasticity than in the convection- 
dominated case. The increasing probability of explosion with increasing heating parameter 
in the SASI-dominated case undermines the paradigmatic notion of a critical luminosity 
for explosion (Burrows & Goshy 1993), or at least modifies it to a smeared-out ‘threshold 
region.]^ 


^The increase of the SASI-dominated luminosity threshold with higher resolution reported by 


Fernandez 
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Is the SASI a potentially important contributor to supernova explosions? With more 
complete physics, it has tended to appear in cases that do not explode, leading some to 
believe that the SASI is a curiosity that does not directly contribute to explosions in nature. 
However, our results on the heating efficiencies of convection-dominated vs. SASI-dominated 
explosions motivate an open mind, particularly given the historical lack of robustness of 
the delayed neutrino-driven mechanism. At first blush, the lower efficiencies of the SASI- 
dominated case in the lower panel of Figure may seem less conducive to explosions, for 
we usually hear about multidimensional effects increasing heating efficiency. That may be 
true of convection. But our intuition about neutrino heating efficiency can be turned on its 
head: if a model can explode with lower efficiency, it requires less from the neutrino heating, 
making the neutrino-driven mechanism more robust]^ 


We have artificially induced the SASI by ‘turning off’ our mock-up of nuclear disso¬ 
ciation]^ but there are also physical possibilities for making the SASI more likely. Several 
actively-studied aspects of progenitor structure affect the ease with which the SASI develops 
(Hanke et ah 2013; Couch & Ott 2015; Muller & Janka 2015). The equation of state at 
low density affects the shock jump conditions, and at high density affects the neutron star 
radius and speed of contraction. Neutrino transport impacts the stalled shock radius and 
the Foglizzo number y; and note for instance that a change in neutrino opacities produced 
stronger SASI activity in a specific case (Melson et al. 2015a). We should remain open to 
such possibilities conducive to the SASI in searching for a robust explosion mechanism, for 
nature may take advantage of its ability to do some of the ‘heavy lifting.’ 


We thank Rodrigo Fernandez for answering questions about his models. We thank Eirik 
Endeve and Anthony Mezzacappa for discussions and ongoing collaboration. This material 
is based upon work supported by the U.S. Department of Energy, Office of Science, Office of 
Advanced Scientific Computing Research and Office of Nuclear Physics. This research used 
resources of the Joint Institute for Computational Sciences at the University of Tennessee. 


(2015), while physically plausible, requires reassessment with more realizations. Such a study would also 
indicate how the magnitude of threshold smearing—roughly 20% in the present study—may depend on 
resolution. 


^See 


Couch & Ott 


(2015) on this ‘inverted’ perspective in connection with turbulence. 


^Another artificial means of investigating explosions—increasing neutrino luminosities beyond nominally 
physical values—tends to suppress the SASI, perhaps contributing to prejudice against it. 
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